This notebook contains the code of the paper "Bayesian Calibration of Imperfect Computer Models using Physics-Informed Priors". The models are fitted in rstan and the code is available in the folder "STAN/Approximations".
knitr::opts_chunk$set(
message=FALSE,
warning = FALSE,
comment = '',
fig.width = 6,
fig.height = 4,
fig.align = 'center'
)
getwd()
[1] "/Users/michais/Desktop/BCPI_codes_new"
# uncomment to install
# install.packages("rstan")
# install.packages("ggplot2")
# install.packages("tidyverse")
library(rstan)
library(ggplot2)
library(tidyverse)
rstan_options(auto_write = TRUE)
options(mc.cores = 3) # allocate 3 cores (for each model we run 3 chains in parallel)
# Numerical simulator of the WK3 model
source("functions/WK2and3_sim_fn.R")
# Load flow data
d = readRDS("Data/Inflow_time.rds")
Rtrue = 1; Ctrue = 1.1; Ztrue = 0.05
flow = d$inflow*0.95
time = d$time
nP = 90 # number of pressure data
nI = 100# number of inflow data
nc = 1 # number of cardiac cycles
nflow = length(flow)
# 1. simulate WK3 data (R=R_2, Z=R_1)
Psim = WK3_simulate(flow = flow, time = time, R = Rtrue, C = Ctrue, Z=Ztrue) # simulate WK3 data for a given flow Q(t)
P_true = Psim
# 2. choose pressure and inflow indices
indP = round(seq(1, nflow, length.out = nP)); indI = round(seq(1, nflow, length.out = nI))
yP_real = Psim[indP]; yI_real = flow[indI] # noise free fimulated pressure and flow
# 3. Add noise
# set.seed(0)
set.seed(1)
Pnoise = rnorm(nP*nc, 0, 4) # sample pressure noise from N(0, 4^2)
Inoise = rnorm(nI*nc, 0, 10) # sample flow noise from N(0,10^2)
yP_real = rep(yP_real,nc)
yI_real = rep(yI_real,nc)
# 4. store individual data in the population matrices
yP = yP_real + Pnoise # add noise
yI = yI_real + Inoise # add noise
tP = time[indP] # corresponding time (synchronized for the two cycles)
tI = time[indI] # corresponding time (synchronized for the two cycles)
data_PI = list(nP=nc*nP, nI=nc*nI, tP=rep(tP,nc), tI=rep(tI,nc), yP=yP, yI=yI, mP=12, mI=10)
WK2_VFE = stan_model('STAN/Approximations/VFE/WK2_delta_VFE.stan')
kp = kmeans(data.frame(x=data_PI$tP), centers = data_PI$mP)
ki = kmeans(data.frame(x=data_PI$tI), centers = data_PI$mI)
init = list("zP" = as.vector(kp$centers), "zI" = as.vector(ki$centers))
op_VFE=optimizing(WK2_VFE, data=data_PI, hessian=FALSE, init = init, verbose=TRUE, seed=0)
Chain 1: Initial log joint probability = -13400
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 19 -962.859 0.0138579 2.3978 1 1 31
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 39 -949.581 0.0272867 2.43968 1 1 56
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 49 -949.485 0.000188747 0.0088334 1 1 68
Chain 1: Optimization terminated normally:
Chain 1: Convergence detected: relative gradient magnitude is below tolerance
op_VFE
$par
rho alpha rho_d alpha_d mu_wk2 sigmaP sigmaI R C zP[1] zP[2] zP[3] zP[4]
0.195718978 4.897341220 0.057918221 0.001593782 97.398040198 9.760864357 63.624564498 0.955552758 1.013489858 0.080072408 0.133867715 0.455149871 0.194511717
zP[5] zP[6] zP[7] zP[8] zP[9] zP[10] zP[11] zP[12] zI[1] zI[2] zI[3] zI[4] zI[5]
0.274922964 0.724878480 0.826984618 0.944283254 0.365017283 0.544772397 0.020043511 0.634880324 0.970363221 0.614334337 0.714632658 0.049797442 0.394214611
zI[6] zI[7] zI[8] zI[9] zI[10]
0.511302773 0.809613912 0.274381891 0.161115333 0.894987597
$value
[1] -949.4846
$return_code
[1] 0
$theta_tilde
rho alpha rho_d alpha_d mu_wk2 sigmaP sigmaI R C zP[1] zP[2] zP[3] zP[4] zP[5] zP[6] zP[7] zP[8]
[1,] 0.195719 4.897341 0.05791822 0.001593782 97.39804 9.760864 63.62456 0.9555528 1.01349 0.08007241 0.1338677 0.4551499 0.1945117 0.274923 0.7248785 0.8269846 0.9442833
zP[9] zP[10] zP[11] zP[12] zI[1] zI[2] zI[3] zI[4] zI[5] zI[6] zI[7] zI[8] zI[9] zI[10]
[1,] 0.3650173 0.5447724 0.02004351 0.6348803 0.9703632 0.6143343 0.7146327 0.04979744 0.3942146 0.5113028 0.8096139 0.2743819 0.1611153 0.8949876
zP_opt_VFE=op_VFE$par[grep("zP",names(op_VFE$par))]
zI_opt_VFE=op_VFE$par[grep("zI",names(op_VFE$par))]
# plot(sort(zP_opt_VFE))
# plot(sort(zI_opt_VFE))
data_PI_Z_VFE= data_PI
data_PI_Z_VFE$zP = zP_opt_VFE
data_PI_Z_VFE$zI = zI_opt_VFE
fit_post_VFE=stan(file='STAN/Approximations/VFE/WK2_delta_VFE_fixed_Z.stan',
data=data_PI_Z_VFE,
chains=3,
iter=1000,
seed=0
)
starting worker pid=67458 on localhost:11354 at 11:30:29.975
starting worker pid=67472 on localhost:11354 at 11:30:30.172
starting worker pid=67486 on localhost:11354 at 11:30:30.365
SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.009392 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 93.92 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.010761 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 107.61 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.009607 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 96.07 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 58.2167 seconds (Warm-up)
Chain 3: 47.0122 seconds (Sampling)
Chain 3: 105.229 seconds (Total)
Chain 3:
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 60.1515 seconds (Warm-up)
Chain 2: 47.3321 seconds (Sampling)
Chain 2: 107.484 seconds (Total)
Chain 2:
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 69.9765 seconds (Warm-up)
Chain 1: 51.3671 seconds (Sampling)
Chain 1: 121.344 seconds (Total)
Chain 1:
# stan_hist(fit_post_VFE)
stan_trace(fit_post_VFE)
zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))
zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))
fit_post_FITC=stan(file='STAN/Approximations/FITC/WK2_delta_FITC_fixed_Z.stan',
data=data_PI_Z_FITC,
chains=3,
iter=1000,
seed=0
)
starting worker pid=67510 on localhost:11354 at 11:32:59.213
starting worker pid=67524 on localhost:11354 at 11:32:59.413
starting worker pid=67538 on localhost:11354 at 11:32:59.614
SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.007932 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 79.32 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.007228 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 72.28 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.006141 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 61.41 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 60.4329 seconds (Warm-up)
Chain 1: 30.205 seconds (Sampling)
Chain 1: 90.6379 seconds (Total)
Chain 1:
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 56.8 seconds (Warm-up)
Chain 3: 38.3528 seconds (Sampling)
Chain 3: 95.1527 seconds (Total)
Chain 3:
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 62.991 seconds (Warm-up)
Chain 2: 38.5151 seconds (Sampling)
Chain 2: 101.506 seconds (Total)
Chain 2:
# stan_hist(fit_post_FITC)
stan_trace(fit_post_FITC)
post_VFE = data.frame(rstan::extract(fit_post_VFE))
post_FITC = data.frame(rstan::extract(fit_post_FITC))
pr = c("R", "C", "sigmaP", "sigmaI")
pVFE = as.vector(as.matrix(post_VFE[,pr]))
pFITC = as.vector(as.matrix(post_FITC[,pr]))
Ns = nrow(post_VFE)
df_plot_post = data.frame(post = c(pVFE, pFITC), par = rep(rep(pr, each = Ns),2), model = c(rep("VFE",length(pr)*Ns), rep("FITC",length(pr)*Ns)))
mod_dat = df_plot_post%>%
mutate(par = recode(par,
"R" = "R",
"C" = "C",
"sigmaP" = "sigma[P]",
"sigmaI" = "sigma[Q]"
))
df_true_par = data.frame(post=c(Rtrue+Ztrue, Ctrue, 4, 10),par = c("R", "C", "sigma[P]", "sigma[Q]"))
set_lim = data.frame(x=c(0.5,3.1),y=c(600, 600))
df_point_est = data.frame(val=c(op_VFE$par[pr],op_FITC$par[pr]),par = rep(df_true_par$par,2), model = rep(rep(c("VFE","FITC"),each=4),2))
pl_post = ggplot()+
geom_histogram(data = mod_dat, aes(x=post, fill = par), color="black",bins = 60)+
facet_grid(model~par, scales = "free", labeller = labeller(par = label_parsed))+
geom_point(data = set_lim, aes(x=x,y=y), color = "white",alpha=0.000001, size=0.000001)+ # set limits on R and C
geom_vline(aes(xintercept = post, linetype = "true"), data=df_true_par, size=0.8)+
geom_vline(aes(xintercept = val, linetype = "map"), data=df_point_est, size=0.8)+
theme_bw()+
theme(#legend.position = "none",
legend.title = element_blank(),
legend.position="bottom",
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
strip.text.x = element_text(size = 13),
strip.text.y = element_text(size = 13))+
xlab("") + ylab("")+
scale_fill_manual(
breaks=c("R", "C", "sigma[P]", "sigma[Q]"),
values=c("#F8766D","#00BE67","white", "white"),guide = "none")
(pl_post=pl_post + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))
ggsave("figures/Appr_delta_post.pdf", plot = pl_post, width = 20, height = 12, units = "cm")
nP_pred = 40
ind_P_pred = round(seq(1,length(time),length.out = nP_pred))
tP_pred = time[ind_P_pred]
data_PI_Z_FITC$nP_pred = nP_pred
data_PI_Z_FITC$tP_pred = tP_pred
data_PI_Z_FITC$nI_pred = nP_pred
data_PI_Z_FITC$tI_pred = tP_pred
N_samples = nrow(post_FITC)
data_post_FITC = list(alpha=post_FITC$alpha, rho=post_FITC$rho, alpha_d=post_FITC$alpha_d
, rho_d=post_FITC$rho_d, sigmaP=post_FITC$sigmaP, sigmaI=post_FITC$sigmaI
, R=post_FITC$R, C=post_FITC$C, N_samples=N_samples
)
data_pred_FITC = c(data_PI_Z_FITC, data_post_FITC)
pred_FITC = stan(file = 'STAN/Approximations/FITC/FITC_delta_predictions.stan',
data = data_pred_FITC,
chains = 1, iter = 1, seed=123,
algorithm = "Fixed_param")
SAMPLING FOR MODEL 'FITC_delta_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0 seconds (Warm-up)
Chain 1: 3.00098 seconds (Sampling)
Chain 1: 3.00098 seconds (Total)
Chain 1:
post_mu_CIs.fn = function(post_pred, ci = c(0.05,0.95), time = tP_pred){
pp = rstan::extract(post_pred)
dfP = pp$y_P[1,,]
Psmr = data.frame(mean = colMeans(dfP),
lower = apply(dfP, 2, quantile, probs = ci[1]),
upper = apply(dfP, 2, quantile, probs = ci[2]),
time=time)
dfI = pp$y_I[1,,]
Ismr = data.frame(mean = colMeans(dfI),
lower = apply(dfI, 2, quantile, probs = ci[1]),
upper = apply(dfI, 2, quantile, probs = ci[2]),
time=time)
return(list(Psmr=Psmr, Ismr=Ismr))
}
data_post_VFE = list(alpha=post_VFE$alpha, rho=post_VFE$rho, alpha_d=post_VFE$alpha_d
, rho_d=post_VFE$rho_d, sigmaP=post_VFE$sigmaP, sigmaI=post_VFE$sigmaI
, R=post_VFE$R, C=post_VFE$C, N_samples=N_samples
)
data_pred_VFE = c(data_PI_Z_VFE, data_post_VFE)
data_pred_VFE$nP_pred = nP_pred
data_pred_VFE$tP_pred = tP_pred
data_pred_VFE$nI_pred = nP_pred
data_pred_VFE$tI_pred = tP_pred
pred_VFE = stan(file = 'STAN/Approximations/VFE/WK2_delta_VFE_predictions.stan',
data = data_pred_VFE,
chains = 1, iter = 1, seed=123,
algorithm = "Fixed_param")
SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0 seconds (Warm-up)
Chain 1: 2.47053 seconds (Sampling)
Chain 1: 2.47053 seconds (Total)
Chain 1:
pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pp_FITC = rbind(post_mu_CIs.fn(post_pred=pred_FITC)$Psmr,post_mu_CIs.fn(post_pred=pred_FITC)$Ismr)
pp_FITC$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_FITC$model = "FITC"
pred_df = rbind(pp_VFE, pp_FITC)
# head(pred_df)
We observe that the VFE model can severely overestimate the noise and therefore result in underfitting. A remedy to this problem is to fix the noise parameter of the functions \(P(t)\) and \(Q(t),\) (\(\sigma_P\) and \(\sigma_I\)). A possible solution for obtaining estimates for the noise parameters is to fit a standard GP model for each dataset \((y_P,t_P)\) and \(y_Q, t_Q\) indepentently and obtain MLE estimates via maximizing the marginal log-likelihood.
df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_FITC = data.frame(z=data_PI_Z_FITC$zP, y = rep(60, data_PI_Z_FITC$mP), model = "FITC", output = "pressure (mmHg)")
df_zI_FITC = data.frame(z=data_PI_Z_FITC$zI, y = rep(-100, data_PI_Z_FITC$mI), model = "FITC", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE, df_zP_FITC, df_zI_FITC)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
facet_grid(output~model,scales = "free")+
geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
scale_fill_manual("",values=c("90% CI" = "grey12"))+
theme_bw()+xlab("time (sec)")+ylab("")+
scale_shape_manual("", values = c("Z" = 4))+
theme(#legend.position = "none",
legend.title = element_blank(),
legend.position="bottom",
strip.text.x = element_text(size = 13),
strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))
ggsave("figures/Appr_delta_pred.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")
nsP = 25
indP = seq(1,data_PI$nP,length.out = nsP)
data_sample_P = list(N=nsP, x = data_PI$tP[indP], y = data_PI$yP[indP])
data_sample_I = list(N=nsP, x = data_PI$tI[indP], y = data_PI$yI[indP])
GP = stan_model('STAN/Approximations/GP_full/GP.stan')
GP_MLE_P=optimizing(GP, data=data_sample_P, hessian=FALSE, verbose=TRUE,seed=0)
Chain 1: Initial log joint probability = -84.3916
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 18 -65.767 0.00556334 7.66548e-05 1 1 30
Chain 1: Optimization terminated normally:
Chain 1: Convergence detected: relative gradient magnitude is below tolerance
GP_MLE_P
$par
rho alpha sigma
0.2147435 76.0190356 3.3560127
$value
[1] -65.767
$return_code
[1] 0
$theta_tilde
rho alpha sigma
[1,] 0.2147435 76.01904 3.356013
sigma_P_MLE = GP_MLE_P$par["sigma"]
GP_MLE_I=optimizing(GP, data=data_sample_I, hessian=FALSE, verbose=TRUE,seed=0)
Chain 1: Initial log joint probability = -656.901
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 19 -111.136 0.558199 0.622198 1 1 26
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 39 -111.126 0.321616 0.000908148 0.4931 0.04931 55
Chain 1: Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
Chain 1: 41 -111.126 0.373032 0.00159772 0.1782 1 58
Chain 1: Optimization terminated normally:
Chain 1: Convergence detected: relative gradient magnitude is below tolerance
GP_MLE_I
$par
rho alpha sigma
0.07883712 99.99999925 11.38052032
$value
[1] -111.1259
$return_code
[1] 0
$theta_tilde
rho alpha sigma
[1,] 0.07883712 100 11.38052
sigma_I_MLE = GP_MLE_I$par["sigma"]
data_pred_VFE_MLE = data_pred_VFE
data_pred_VFE_MLE$sigmaP = sigma_P_MLE
data_pred_VFE_MLE$sigmaI = sigma_I_MLE
pred_VFE_MLE = stan(file = 'STAN/Approximations/VFE/MLE_sigma/WK2_delta_VFE_predictions.stan',
data = data_pred_VFE_MLE ,
chains = 1, iter = 1, seed=123,
algorithm = "Fixed_param")
SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0 seconds (Warm-up)
Chain 1: 2.4886 seconds (Sampling)
Chain 1: 2.4886 seconds (Total)
Chain 1:
pp_VFE_MLE = rbind(post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Ismr)
pp_VFE_MLE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE_MLE$model = "VFE fixed noise"
pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pred_df = rbind(pp_VFE_MLE,pp_VFE)
df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zP, y = rep(60, data_pred_VFE_MLE$mP), model = "VFE fixed noise", output = "pressure (mmHg)")
df_zI_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zI, y = rep(-100, data_pred_VFE_MLE$mI), model = "VFE fixed noise", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE,df_zP_VFE_MLE,df_zI_VFE_MLE)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
facet_grid(output~model,scales = "free")+
geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
scale_fill_manual("",values=c("90% CI" = "grey12"))+
theme_bw()+xlab("time (sec)")+ylab("")+
scale_shape_manual("", values = c("Z" = 4))+
theme(#legend.position = "none",
legend.title = element_blank(),
legend.position="bottom",
strip.text.x = element_text(size = 13),
strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))
ggsave("figures/Appr_delta_pred_noise.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")